3.Image Analysis
library(data.table)
library(anytime)
library(stringr)
library(dplyr)
library(gridExtra)
library(ggplot2)
library(plotly)
library(jpeg)
library(fields)
library(formatR)
img = readJPEG("/Users/ilkerkurtulus/Downloads/IMG_7955.jpg")
str(img)
## num [1:512, 1:512, 1:3] 0.424 0.4 0.388 0.404 0.4 ...
dim(img)
## [1] 512 512 3
if (exists("rasterImage")) {
plot(1:2, type = "n")
rasterImage(img, 1, 1, 2, 2)
}
if (exists("rasterImage")) {
plot(1:4, type = "n")
rasterImage(img[, , 1], 1, 1, 2, 2)
rasterImage(img[, , 2], 2, 1, 3, 2)
rasterImage(img[, , 3], 3, 1, 4, 2)
}
noise = runif(512 * 512 * 3, min = 0, max = 1)
dim(noise) = c(512, 512, 3)
nimg = noise + img
nimg_scaled = (nimg - min(nimg))/(max(nimg) - min(nimg))
if (exists("rasterImage")) {
plot(1:2, type = "n")
rasterImage(nimg_scaled, 1, 1, 2, 2)
}
if (exists("rasterImage")) {
plot(1:4, type = "n")
rasterImage(nimg_scaled[, , 1], 1, 1, 2, 2)
rasterImage(nimg_scaled[, , 2], 2, 1, 3, 2)
rasterImage(nimg_scaled[, , 3], 3, 1, 4, 2)
}
to_greyscale = function(img) {
img[, , 1] * 0.21 + img[, , 2] * 0.72 + img[, , 3] * 0.07
}
grey_nimg = to_greyscale(nimg_scaled)
pca_img = prcomp(grey_nimg, center = TRUE, scale. = TRUE)
eigs = pca_img$sdev^2
exp_var_ratio = eigs/sum(eigs)
cum_exp_var_ratio = cumsum(exp_var_ratio)
plot(cum_exp_var_ratio)
Even the image has 512 dimensions we can explain it well with 250 components with 0.95 confidence interval. However its still very large dimension for statistical inference. b) Construct patch matrix:
extract_path = function(data, i, j, n) {
dw = (n - 1)/2
as.vector(t(data[(i - dw):(i + dw), (j - dw):(j + dw)]))
}
tmp = c()
for (i in 2:511) {
for (j in 2:511) {
x = extract_path(grey_nimg, i, j, 3)
tmp = c(tmp, x)
}
}
res = matrix(tmp, nrow = 260100, byrow = TRUE)
Apply pca and plot:
pca_res = prcomp(res, center = TRUE, scale. = TRUE)
# 1st component
pca_res1d = pca_res$x[, 1]
pca_mat = matrix(pca_res1d, nrow = 510, ncol = 510, byrow = TRUE)
pca_mat_scaled = (pca_mat - min(pca_mat))/(max(pca_mat) - min(pca_mat))
if (exists("rasterImage")) {
plot(1:2, type = "n")
rasterImage(pca_mat_scaled, 1, 1, 2, 2)
}
# 2nd component
pca_res1d = pca_res$x[, 2]
pca_mat = matrix(pca_res1d, nrow = 510, ncol = 510, byrow = TRUE)
pca_mat_scaled = (pca_mat - min(pca_mat))/(max(pca_mat) - min(pca_mat))
if (exists("rasterImage")) {
plot(1:2, type = "n")
rasterImage(pca_mat_scaled, 1, 1, 2, 2)
}
# 3rd component
pca_res1d = pca_res$x[, 3]
pca_mat = matrix(pca_res1d, nrow = 510, ncol = 510, byrow = TRUE)
pca_mat_scaled = (pca_mat - min(pca_mat))/(max(pca_mat) - min(pca_mat))
if (exists("rasterImage")) {
plot(1:2, type = "n")
rasterImage(pca_mat_scaled, 1, 1, 2, 2)
}
ev1 = pca_res$rotation[, 1]
ev_scaled1 = (ev1 - min(ev1))/(max(ev1) - min(ev1))
mat_ev1 = matrix(ev_scaled1, nrow = 3, ncol = 3, byrow = TRUE)
ev2 = pca_res$rotation[, 2]
ev_scaled2 = (ev2 - min(ev2))/(max(ev2) - min(ev2))
mat_ev2 = matrix(ev_scaled2, nrow = 3, ncol = 3, byrow = TRUE)
ev3 = pca_res$rotation[, 3]
ev_scaled3 = (ev3 - min(ev3))/(max(ev3) - min(ev3))
mat_ev3 = matrix(ev_scaled3, nrow = 3, ncol = 3, byrow = TRUE)
if (exists("rasterImage")) {
plot(1:4, type = "n")
rasterImage(mat_ev1, 1, 1, 2, 2)
rasterImage(mat_ev2, 2, 1, 3, 2)
rasterImage(mat_ev3, 3, 1, 4, 2)
}